
dir.create("PSID_eventstudy")

rawdata<-as.data.table(read.dta13("psid_cleaned_headanysp.dta"))
rawdata[,onset_year := max((age==onset_age)*year),by=newid]
rawdata[,anyonset_year := max((age==anyonset_age)*year),by=newid]
rawdata[is.na(onset_year) | onset_year ==0,onset_year := Inf]
rawdata[is.na(anyonset_year) | anyonset_year ==0,anyonset_year := Inf]
rawdata[is.na(anyonset_age),anyonset_age:=Inf]
rawdata[is.na(onset_age),onset_age:=Inf]
rawdata[,pre96:= as.numeric(year<1996)]

rawdata[,agebin:=floor(age/10)*10]

rawdata[year<1996,baseyear:=-1]
rawdata[year>=1996,baseyear:=-2]

setorder(rawdata,newid,year)
rawdata[,empyet:=cumsum(emp)>=1,by=newid]
rawdata[,empyet:=empyet == 1 & sample_analytic == 1,by=newid]
rawdata<-rawdata[!is.na(onset_age),]
rawdata[,N:=1]
rawdata[,time:=age-onset_age]

rawdata[,empyet1:=max(empyet == 1 & (time == 0)), by=newid]
rawdata[,empyet2:=max(empyet == 1 & (time == baseyear)), by=newid]

rawdata[,married1:=max(marnow == 1 & (time == 0)), by = newid]
rawdata[,married1b:=max(marit == 1 & (time == 0)), by = newid]
rawdata[,married2:=max(marnow == 1 & (time == baseyear)), by = newid]
rawdata[,married2b:=max(marit == 1 & (time == baseyear)), by = newid]
rawdata[married1b==0&married2b==0,mardynamic:=0]
rawdata[married1b==1&married2b==0,mardynamic:=1] #gets married
rawdata[married1b==0&married2b==1,mardynamic:=2] #gets divorced
rawdata[married1b==1&married2b==1,mardynamic:=3] #always married

feols(anyDI ~ married1, data = rawdata[time>= 0 & time <= 4 & empyet1==1,],cluster="newid")
feols(anyDI ~ married1b, data = rawdata[time>= 0 & time <= 4 & empyet1==1,],cluster="newid")
feols(anyDI ~ married2, data = rawdata[time>= 0 & time <= 4 & empyet2==1,],cluster="newid")
feols(anyDI ~ married2b, data = rawdata[time>= 0 & time <= 4 & empyet2==1,],cluster="newid")

rawdata[,married:=married1b]
#rawdata[,empyet:=empyet1]

rawdata[,ownhome:=max(home_equity >1 & (time <= baseyear),na.rm=TRUE),by=newid]
rawdata[,highedu:=educ>=3]
rawdata[,ageold:=onset_age>=45]
for(agelev in 0:1){
  spwagemed<-rawdata[ageold==agelev & sp_wage>0&!is.na(sp_wage)&time<=0&married==1,quantile(sp_wage,0.5)]
  rawdata[ageold==agelev&sp_wage>0&!is.na(sp_wage)&time<=0&married==1,highspwage:=as.numeric(sp_wage >= spwagemed)]
}
rawdata[,highspwage:=max(highspwage,na.rm=TRUE),by=newid]
rawdata[is.infinite(highspwage),highspwage:=NA]
rawdata[married==0,highspwage:=NA]
rawdata[married==1,workingspouse:=1]
rawdata[is.na(highspwage)&married==1,workingspouse:=0]

rawdata[,all:=factor(1,levels=c(1,2))]


rawdata[,marhome:=married==1 &ownhome==1]
rawdata[married!=ownhome,marhome:=NA]

rawdata[year<=2001-1 & occupation%in%1:195, occ_group := "Professionals"]
rawdata[year<=2001-1 & occupation%in%201:245, occ_group := "Managers"]
rawdata[year<=2001-1 & occupation%in%260:285, occ_group := "Sales"]
rawdata[year<=2001-1 & occupation%in%301:395, occ_group := "Clerical"]
rawdata[year<=2001-1 & occupation%in%c(401:600,601:695), occ_group := "Craftsmen"]# and Trades"]
rawdata[year<=2001-1 & occupation%in%701:715, occ_group := "Transportation"]
rawdata[year<=2001-1 & occupation%in%740:824, occ_group := "Labor"]# and Farm"]
rawdata[year<=2001-1 & occupation%in%901:984, occ_group := "Service"]# and Household"]

rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%1:43, occ_group := "Managers"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%50:395, occ_group := "Professionals"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%c(400:465), occ_group := "Service"]# and Household"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%470:496, occ_group := "Sales"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%500:593, occ_group := "Clerical"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%c(600:613,770:896), occ_group := "Labor"]# and Farm"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%c(620:676, 700:762, 680:694), occ_group := "Craftsmen"]# and Trades"]
rawdata[year>2001-1 & year <= 2015 - 1 & occupation%in%900:975, occ_group := "Transportation"]


rawdata[year>=2017-1 & occupation%in%10:430, occ_group := "Managers"]
rawdata[year>=2017-1 & occupation%in%500:3955, occ_group := "Professionals"]
rawdata[year>=2017-1 & occupation%in%c(4000:4650), occ_group := "Service"]# and Household"]
rawdata[year>=2017-1 & occupation%in%4700:4965, occ_group := "Sales"]
rawdata[year>=2017-1 & occupation%in%5000:5940, occ_group := "Clerical"]
rawdata[year>=2017-1 & occupation%in%c(6000:6130,7700:8965), occ_group := "Labor"]# and Farm"]
rawdata[year>=2017-1 & occupation%in%c(6200:6940, 7000:7630), occ_group := "Craftsmen"]# and Trades"]
rawdata[year>=2017-1 & occupation%in%9000:9750, occ_group := "Transportation"]

setorder(rawdata,newid,year)
rawdata[,bluecollar:=occ_group%in%c("Service","Labor","Craftsmen","Transportation")]
rawdata[DS!=0,bluecollar:=NA]

for(occ in unique(rawdata[!is.na(occ_group),]$occ_group)) {
  rawdata[,paste0("occ_group_",eval(occ)):= nafill(as.numeric(occ_group==occ),type="locf"),by=newid]
}

wealthcutoff<-rawdata[year==onset_year+baseyear & year >= 1998&!is.na(w_no_h),quantile(w_no_h,probs=0.75)]

rawdata[year>=1998,rich:=w_no_h>=wealthcutoff]
rawdata[year==onset_year+baseyear & year >= 1998,prerich:=rich]
rawdata[,temp:=all(is.na(prerich)),by=newid]
rawdata[,prerich:=max(prerich,na.rm=TRUE),by=newid]
rawdata[temp==TRUE,prerich:=NA]
rawdata[,temp:=NULL]

#carry forward last non-missing obs:
rawdata[,bluecollar:=nafill(as.numeric(bluecollar),type="locf"),by=newid]
rawdata[,N:=1]
rawdata[,n:=1]
rawdata[,health:=DS+1]

if(eventmeans==TRUE){
  rawdata[,empyet:=empyet2]
  rawdata_mod<-rawdata#[onset_age_hsev>onset_age + 5 | is.na(onset_age_hsev),]
  rawdata<-rawdata[onset_age_hsev<onset_age + 5,]
  

}
